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In this paper, we show that using measurements for different frequencies, and using ultrasound 
localized perturbations it is possible to extend the method of the imaging by elastic deforma- 
tion developed by Ammari and al. [Electrical Impedance Tomography by Elastic Deformation 
SIAM J. Appl. Math. , 68(6), (2008), 1557-1573.] to problems of the form 

div(7V«) + k^qu =0 in 57, 
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and to reconstruct by a perturbation method both 7 and q, provided that 7 is coercive and 
k is not a resonant frequency. 
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1. Introduction and Notations 

In the recent years, a lot of attention has been devoted to the reconstruction of 
physical parameters of partial differential equations from electromagnetic mea- 
surements. In the case of electrical impedance tomography (EIT) it is well known 
that the detection of the conductivity from boundary measurements is a very ill- 
conditioned problem. This drawback has limited its use so far to anomaly detection. 
In a recent work, Ammari et al. have shown that combining these measurements 
with simultaneous localized ultrasonic perturbations allows to recover the con- 
ductivity with great precision. The purpose of this work is to show that such an 
approach can be generalized successfully to the study of Helmholtz type problems. 
In what follows we use the following notations: 

• is a smooth domain in IR" with a regular boundary denoted by 90, 

• X is a point in il, 

• Vt' = {x &Vt\ dist(x,(90) > do > 0} represents the interior points of 0, 

• w C 0' is the region of the localization of the ultrasound perturbations, which 
is supposed to be small compared to the size of 0', 

• \w\ is the volume of w, 
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• 1^ denotes the characteristic function corresponding to the set w, i.e., the func- 
tion which takes the value 1 on the set and the value outside, 

• z £ w is the centre point of the region of the ultrasound perturbation, 

• fcj > is a frequency, 

• 7(x) is the conductivity and is a scalar real- valued function such 
that < Co < 7(x) < Co for all x G 

• q{x) is the permittivity and is a scalar real-valued function such 
that < Co < q{x) < Cq for all x € O, 

• u{x) is the potential induced on the boundary by the electromagnetic field (p 
in the absence of ultrasonic perturbations (n(x) and (p{x) are complex-valued 
functions) , 

• Uw is the perturbed potential field induced on the boundary by the electromag- 
netic field (f in the presence of ultrasonic perturbations localized in the domain 
w {uw is a complex- valued function), 

• A is the amplitude of the ultrasonic perturbation, 

• 'yw{x) is the perturbed conductivity (real- valued positive bounded function), 

• 7 is the value of the perturbed conductivity 7^, in the area w of the perturbation 
(real- valued positive bounded function), 

• qw{x) is the perturbed permittivity (real- valued positive bounded function), 

• q is the value of the perturbed permittivity in the area w of the perturbation 
(real- valued positive bounded function), 

• and m^„ are the polarization tensors, 

• N^^q{x,z) is the Neumann function for the operator div(7(2;)V2:) + q{x) in VL 
corresponding to a Dirac mass at 2, 

• W^{^Y) is the Sobolev space of the functions u{x) such that u € Lf^{VL) and 
Vu € L^{VL), 

• for the complex-valued function tt, the function u denotes its complex- 
conjugated. 

The problem we consider is the following. Let 7 € C^(il) and q € C^{^) be 
bounded scalar real- valued functions (see the list of notations). For i = 1,2, let 
Ui € H'^{^) be such that 

div(7Vnj) + kfqui =0 in J^, (1) 
du ' 

7—^ = on (90. (2) 

Of 

The well-posedness of this problem requires that kf is not an eigenvalue of the 
generalized eigenvalue problem 

— div(7Vn) = \qu in il, (3) 

7--— =0 on dVL. 
ov 

It is well known that this problem admits a countable number of eigenmodes, 
with no accumulation point, and that each eigenvalue as a finite multiplicity. We 
will assume that ki and /c2 do not correspond to eigenvalues of problem ([3]). The 
generalization of the method introduced in [if is the following. For frequency fcj 
being fixed, we measure the potential Uj, solution of problem ([I])-(l2|), on 317. 

Assume now that ultrasonic waves are localized around a point z € fi, creating a 
local change in the physical parameters of the medium. Further, we suppose that 
q and 7 are known close to the boundary of the domain, so that ultrasonic probing 
is limited to interior points x in Vt' (see the list of notations), where do is very large 
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compared to the radius of the spot of the ultrasonic perturbation. 

We suppose that this deformation affects 7 and q linearly with respect to the 
amphtude of the ultrasonic signal. Such an assumption is reasonable if the am- 
plitude is not too large. Thus, when the electric potential is measured while the 
ultrasonic perturbation is enforced, the equation for the potential is 



div(7^„Vni,^) + kfqyjUi^u, = in $7, (4) 
dv 



7-^^ = '^i on dVt. (5) 



with 



7t« =7 + lt«(7A-7), (6) 
qw = q + lyo{q\- q), (7) 

where A is the amplitude of the ultrasonic perturbation given by the ratio of the 
perturbed volume Vw oi w over the unperturbed one Vw (see [H). In other words 



ln,KX) - I x{x)^{x), xGw ^'"^^^ ~ \ Xix)q{x), xGw 

where A(x) = VS/V^ is a known function. 

The analysis of the change of the Neumann-to-Dirichlet map as a result of elec- 
tromagnetic perturbation of small volume follows [l|l . The main differences between 
the case of the conductivity equation considered in [l|] and our case of the Helmholtz 
equation are the following: this time the boundary data (p and the solutions Ui are 
complex- valued functions in our case while they are real in [![) and in our case 
we need to reconstruct simultaneous two coupled real-valued parameters 7 and q. 
Therefore we expand the main ideas of [ll] to our case (see Section [2|) . The choice 
of real 7 and q implies the existence of eigenfrequencies (see problem ([3])) and this 
gives an additional difficulty in numeric reconstruction. The case of complex 7 and 
q which allows to avoid the resonances, will be considered in 

The signature of the perturbations on boundary measurements can be measured 
by the change of energy on the boundary, namely 



{uw—u)(pda = \w\ 
an 



7A 



( ) (7A - 7)Vu(z) • Vu{z) - k^{q\ - q)u{z) ■ u{z) 



(8) 

Assuming the perturbed region is 9, bo-ll, tli6 polcirizcLtioii tensor ( ~r~ J is 3, 
sccilctr, 

= (7A-7)/(7A + 7). 

Therefore, for a localized perturbation focused at a point z, we read the following 
data (rescaled by the volume) 

D,{\) = 7|Vn(z)p ^y ~^^ - k\\u{z)\\^\ - 1). (9) 



We notice that the data Dz{X) from Q can be measured for given A and k thanks 
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to the identity: 

= 7-7/ i'^w - u)Tpda. 
\w\ Jan 

The parameters ^(z) and ^(2;) are unknown, but the amplitude A is known. 
Varying the position of locahzation, we are able to recover this localized internal 
data everywhere inside the domain. Thanks to the following lemma 

Lemma 1.1: If the data is known for four distinct values of X, chosen inde- 
pendently of ^ and q, then one can recover ^{z)\Vu{z)\'^ and q{z)\u{z)\'^ . 

we can find directly the functions J{z) = j{z)\Vu{z)\'^ and j{z) = q(z)\u{z)\'^ for 
the unique solution u of problem ([I])- ([2]). 

The proof of Lemma ll. II is simply a study of functions of one variable, which is 
detailed in Appendix lAl 

The rest of the paper is organized as follows: in Section [2] we prove formula (j8]), 
in SectionOwe describe a reconstruction method by perturbations and in Section[4] 
we give and analyse our numeric results, obtained for two different frequencies and 
one boundary data in the form of a plane wave. 



2. Proof of asymptotic expansion ([8]) 

We suppose that k"^ do not correspond to eigenvalues of problem ([3]). To prove the 
asymptotic expansion ([SD, we first need the following Proposition: 

Proposition 2.1: We have the following identities 
div(7^V(u^ -u)) + k'^qwiuw - u) 

= - div(l^(7^ - 7)Vu) - A:^l^(g^ - q)u, (11) 
div(7V(u^ - u)) + k^q{uu, - u) 

= - div(l^(7^„ - '-f)Vuu,) - k^lw{qw - q)uw, (12) 

Thanks to Proposition [2T1 we can estimate the difference between the perturbed 
and unperturbed solutions Uw—u'va L2{^) by a norm of u in the perturbed region 
w and by a power of the small volume \w\ bigger than 0.5. 

Lemma 2.2: Suppose that C i?" contains a subset of^l' C fl of class C^, such 
that dist(0',9r2) > do > 0, and such that w C f^'. Let q, 'j Loo{^) be positive 
functions, satisfying < cq < q{x),^{x) < Cq < +00 a. e. x € Q, and k"^ is not 
a Neumann eigenvalue for problem (0^. Then for the functions and u € W^^ 
verifying Eq. [10]) and Eq. ( fiTj) we have 

\\uw - uWn^iyi) < C|u;|2 |n|tyj^(^). (13) 

Therefore, thanks to relation for m = max{2, n} and all k satisfying < k < 
^ there exists a positive constant C > depending only on , do, cq, and Cq, 
such that 



\\u. 



•w 



-u\\u{n)<C\w\''''"''Hw^{w)- (14) 
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Proof. The proof of estimate (I13p follows the proofs of Lemma 15.1 and Proposi- 
tion 15.2 from [2]. Indeed, as soon as q{x)/^{x)k'^ is not an eigenvalue for the opera- 
tor —A in L2{^) with the homogeneous Neumann boundary condition, in our case 
problem ([I])-([2]) has a unique weak solution u in H^{Q) (for every (p G H~^^'^{dQ)). 
As is uniformly continuous and uniformly coercive on x H^. The embedding 
H^{Q) <S: L2{^) is still compact because dQ € and 17 is compact. 

For passing to the perturbed problem, we change 6 on w (repeat the procedure 
from [4]) and obtain with the help of relation (llip the desired estimate (I13p . 

Let us prove estimate ([Till . Select v as the solution to 



div(7Vt') -|- k'^qv = Uu, — u 
^^1 n 



For this v we have Hi^Hz/a^Q) < C\\uw — "ullLaCn)' ^'^d 



m,„ — u\ 



< 



< 



dx = — 'j{x)V{uuj — u)\/vdx + / uP'q{x){uu, — u)vdx 
n Jn Jfi 

V [d\v{'y{x)'V{uw — u)) + k'^q{x){uw — u)\ dx 

i 

/ V [div(l^(7u_, - ^)Vuw) + k'^l^iqw - q)u^] dx 
Jn 

[div(l^(7^ - 7)Vnu_,) + k^l^iQw - q)uw\ dx 
Jn Jn 



< C 



\SJun,,\'^dx 



Vv\Pdx + { luj^'dx 



Vuwl'^dx] \\v\\h2(^^) + 



in^l^dx ) \\v\\Hi{n) 



Vunjl'^dx] +1 \uy,\'^dx\ \ \\uui - u\\l^(^ci) 



(15) 



provided p, p and q, q are related by | -|- | = 1 and i -|- i = 1. We use Sobolev's 
Embedding Theorem to provide the inclusions C Wp and C Lp. We require 
that q, q > so that I < p, p < For any 1 < q < 2 (see [3, p. 164]) we 

have 



UwWheiiw) < \\Uw - u\\Le,{w) + IfIIl,~(«>) < 

< ( / Idx II u„, — til 



lisW + l^^l'|pllLoe(t") 



< Cl-^l^ ||u||i^(^). 



(16) 
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and for any 1 < q < 2 we obtain 

||Vuu,|1l,{^) < C'glwl'! ||Vu||l^. (17) 
A combination of estimations (|15p . (jl7p and (jl6p yields 

\\uw - u\\L2{n) < C\w\ - + |Vn|i^(^)) 

for any < gr < 2. In other words, ^ < | < = 5 + ^ with = max{2, n}, 
from where we can take q = - = ^ + k for < k < ^. □ 

In addition of estimates ([T3|) and , let us show that the difference u — u.^ can 
be totally described by an integral expression over w. 

Proposition 2.3: Suppose that k"^ is not the Neumann eigenvalue for 
div{'y{x)Vx) + q{x) onw. LetN^q{x,z) be the Neumann function fordiv{'y{x)S/x) + 
q{x) in corresponding to a Dirac mass at z. That is N^g is the solution to 

dW{'y{x)V xN^q{x, z)) + k'^q{x)N^q{x, z) = -5^, in Vt, 
7-3^ =0 on ail. 

Then, by definition of N^q (which is a real function!), the function U defined by 



U{x)= I N^q{x, z)ip{z)da{z) 
Jan 

is the solution of system Therefore, the solutions u and Uw of systems HJ)- 

(0) and gP-(E3j satisfy 

u-Uw){z) = I {'~iw-'~i){z)Vuyj{z)\/ zN^q{z,x)dz+ i k^{q-qw){z)uu,{z)N^q{z,x)dz 



w 



(19) 



Proof. Note that the Neumann function N^q(x, z) is defined as a function of a; G 
Q, for each fixed z ^ Q. Since k"^ is not the Neumann eigenvalue for div(7(x)Va;) + 
q{x) on w, the direct problem ([1]) admits a unique solution u (see [1]). Thus, the 
solution u is represented by the formula 

u{x) = / N^q{x, z)(p{z)da{z). 
Jan 

We notice that 

div{-f{x)Vuy,) + k'^q{x)uw = - div(li„(7^ -7)(x)Vu^) - A;^li„(g'^ - g')(x)n^, (20) 
We multiply relation ([2n|) by A^^^ and integrate over fi: 



(p{z)N^q{x, z)da{z) — / 'y{z)Vuw{z)'VN-yq{x, z)dz + / k q{z)uw{z)N^q{z,x)dz 
an Jq Jn 

/ lw{z){jy, - j){z)\/uy,{z)\/N^q{x,z)dz - / k'^lw{z){qy, - q){z)uy,{z)N^q{x,z)c 
In Jn 
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Therefore, using fQ^Uj^,{z)^{z)V Ng{x, z)da{z) = 0, from the fohowing equahty 

u{z) + / Uwiz) {div{'y{z)VN.yg{x, z)) + k'^q{z)N^g{z,x)) dz = 
Jn 

'i-wiz)ijw - '~l){z)Vuw{z)VN^q{x,z)dz + / k^lw{z){q - qw){z)uw{z)N^q{x, z) 
n Ja 

we obtain Eq. 

Multiplying Eq. ()19p by ^(x) and integrating over dVl, we find 

/ [u - u^)Tpda{x) = / (7u. - 7)(z)Vn^(2;)V^ ( / N^q{z,x)-ip{x)da{x)] dz 
Jan Jw \Jdn / 



+ k {q - qy,){z)uy,{z) [ / N.,g{z,x)ip{x)da{x) \ dz, 
Jw \Jdn / 

which gives 

(« - uj^d.(.) .)(.)V„„(.)V«. + / k% - 

(21) 

Remark 1 : 0] Consider a sequence of sets CC f^. Since the family of functions 
j:^^w^ is bounded in Li{Q), it follows from a combination of the Banach-Alaoglu 
Theorem and Riesz Representation Theorem that we may find a regular, positive 
Borel measure fi, and a subsequence w^^^, with [we^l — )■ 0, such that — )• d/i- 



Finally, thanks to the a priori estimations (|13p . (jl4p and the representation for- 
mula (I21D. we establish the main result: 



Lemma 2.4: Assume that u € Consider a sequence of sets w CC Q 

such that j^^^w converge 
\w\ tends to zero. Then, 



such that j^^^w converges in the sense of measures to a probability measure dfi as 



/ {uu,-u)Tpda = I Mw{^\- -f)\Vu\^dx-k'^\w\ / mw{q\ - q)\u\^dx+0{\w\^^'^). 

J 80. J w J w 

(22) 

The exponent k only depends on Qi, sup^ \qw\, sup^ \'^w\, inff^ \qw\ and inff^ \lw\- 
The remainder term has the form 

\0{\w\^+'')\ < C\w\'-^^\\u\\wM\\'^'t'\\L^m^ 

where C depends only on fii, sup^ \qw\, sup^ |7^|, inf^ \q,u,\ and inf^ |7«,|. Finally, 
with a hypothesis that w is a ball, the polarization tensors and become the 
scalar functions and niw, which are given by 

1 /iA(x)-l\ 1 
Mw = -r—rlw{x) and = -—lyj{x). 

\w\ \^iA(x) + iy \M 

Proof. Suppose that k"^ is not a Neumann eigenvalue for problem ([3|). We have 
relation (j2ip . We are looking for an approximation of the terms of Eq. (j2ip depend- 
ing on Uw by a function depending on u. In the same way as in [l|, we introduce 
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the solution (^^ of the fohowing problem 

J div(7^(2;)VCt„) + k'^qUx)Cw = div(7(2;)V^C) + k^q{x)C in n, 

Corresponding to Cw, we define in the unperturbed case C = x + C+i{x + C), where 
C and C are constants in R"^ for x G R"^. This time all functions Cto, u and Uw 
are complex. The choice of C (C) will be discussed later. Thanks to Lemma 12.21 
for Cu; ~ C we still have an analogue version of Proposition 3.1 of [l|, p. 6]: 

Proposition 2.5: Consider a sequence of sets w CC f^, such that j^^w converges 
in the sense of measures to a probability measure dfi as \w\ tends to zero. Then, 
the corrector converges in the sense of measures to Mjdfj, (Mj is a scalar 

function). Furthermore, it satisfies 

l|V(C«, - C)llL.(n) < C\w\'^ and IIC^ - C||L,(n) < 

where the constants k > and C > depend only on Qi, sup^ \qw\, sup^ {jwl, 
infn and inf^ |7^|. 



The rest of the proof follows the analogous one given in details in This time 
the remaining term is bounded by 



We also remark (see [1[ for the notations) that the choice of ipi = -^u-kri (where 

T] is the standard mollifier) determine the constants C = (Ci, . . . ,Cd) and C = 
(Ci, . . . , Cd) in the definition of the function C,{z): 



Cj + iCj = jT-^ tor a zq € 



which ensures that C,{z)'il) u{z). 
Finally, we deduce 



(n — Uw)fda{x) 
dn 

\w\ [ M^{^^-j){z)\Vu{z)\^dz+\w\ [ k^h!L(q - q^){z)\u{z)\^dz + 0{\w\'+^), 
J m .III} 



^A(a:)+i ) ^ sphere. This proves relation dH]) and 

provide the existence of a known function Dz{X) from □ 
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3. Reconstruction 7 and q by a perturbative method. Numeric algorithm 

We consider the system of Helmholtz equations with different frequencies fci ^ A;2 : 

div(7(x)VufcJ + klq{x)uk, =0 in (23) 
div(7(x)VufcJ + /cig(rE)ufc, = in fi, (24) 
""fci =Uk2 =ip on 9r2. (25) 

The data ip is the Dirichlet data measured as a response to the current ip in absence 
of elastic deformation. We take iJj = Q^^'^'^^'^y/^ ^ which represents a plane wave. 

We use the following formulas 7(x)|Vtifejp = Jk^ix) and ^(x)!^^^^ = jk2{x). 
Thus, we can approximate our problem by system (pSj) and ([271) 

div (^-M^^uk}j + klq{x)uk, =0 in 0, (26) 

div(7(x)VnfcJ + kl^-h^uk, =0 in f^, (27) 

where it is supposed that 

|VnfciP>0 and jn^J^ > for aU x G 0. 

Let us explain the steps of the numeric algorithm. The method uses two sub- 
algorithms to reconstruct 7 for a fixed q (constant for the ultrasound perturbation) 
and to reconstruct q for a fixed 7 (constant). 

First we notice that we have two frequencies ki and k2. 

Step 0. We construct the functions Jk^ and jk^- 

Step 1. We take an initial guess qo and 70. 

Step 2. In the aim of updating first 70 we solve the linear system for chosen q^ 
and 7o and the frequency ki: 

f div(7oV-UfeJ + klqoUk, = 
1 Uk^dn = 'tp 

We obtain the solution of this system which we denote by UQk^- Knowing the 
approximate solution uofc^, we calculate the error on 7: 

Step 3. We verify the condition jii^ofci I < ^precision for a given positive constant 
Cprecisiom which gives the desired order of the precision of the final result. If I-Eofej 
is smaller than Cprecisiom 

we take 7 = 70 and go to Step 5 for the reconstruction of 
g, otherwise we go to Step 4. 

Step 4. We apply the algorithm described in details in Subsection 13.11 to deter- 
mine the correctors and 5uik^ for a fixed qo and to update 70 using formula (|35]) . 

Step 5. In the aim of updating we solve the following linear system with the 
frequency k2 for a chosen go and 70 updated on Step 4: 

f div(7oVufcJ + klqoUk^ = 
yuk^an = ip 
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We obtain the solution of this system which we denote by uq^^. Knowing the 
approximate solution iiofea, we calculate the error on q: 

jk2 

Step 6. We verify the condition leofcj | < eprecision- If \^0k2 1 is smaller than eprecisioiD 
we take q = q^ and finish the algorithm, otherwise we do Step 7. 

Step 7. We apply the algorithm described in details in Subsection l3.2l to determi- 
nate the correctors 5qi and 5uij^^ for a fixed 70 and to update qo using formula (|l0|) . 
Next we go to Step 2. 



3.1. Algorithm of reconstruction of ^ for a constant q 

Step 1. We start from an initial guess 70, and solve the corresponding Dirichlet 
problem for the Helmholtz equation 

div(7o(x)Vuo) + k'^quo = 0, 
uoldn = V'- 

Solving the direct problem for ip = g« Cretan y/x^ obtain uq. 

Step 2. We have seen that our inverse problem is asymptotically approached by 
the direct problem 



div(^^V'u) + uj^qu = 0, in n, 
u = tp, on d^l. 



(28) 



We compute the difference 



and verify 

l-^'ol < C'prec, (30) 



where Cprec is our wished order of the precision. If condition (j30p holds, we finish 
our algorithm and set 7 = 7o- Otherwise we go to the next step. 
Step 3. We use now the expression 

(7o + <57i)|V(no + 5ui)|2 = J{x), 

having the goal to approximate the known function J{x) with the help of the small 
correctors bu\ and We suppose that (5 <C 1 and that 5 max 1 71 1 and 5 max I ui I 

X X 

are of the order of 5. 
By expanding the expression, we obtain 



(1 + 25 



X (1 , /V(Reno)V(Reui) + V(Imtxo)V(Imui)\ ^alV^ipX J{x) 



^^7o(V(Reuo)V(Reui) + V(Imuo)V(Imni)) ^2 |Vni 
~2() 7^ — 7o' 



iVnoP '^iVnol 



2 • 
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We consider only terms of order not smaller than 6: 

. J{x) _7o(V(Reno)V(Re^ii) + V(Im^io)V(Imtxi)) 

= ^ - "° - ■ 

To find the corrector ui = 6ui, we expand the following equation 



11 



div((7o + (^7i)V(tio + dui)) + k^q{uQ + 5ui 



0. 



By considering the terms of order not smaller than 5 and by replacing by 
the approximated formula, we can find ui as the solution of the following problem 



div 



70 V-ui - 2 



Vno 



+ div(£'oV{ti) + div(£^oVuo) + k^qui = 
Ul\dQ. = 0. 



(VReuoVReui + VlmuoVImui^ 



(31) 
(32) 



Let us define 
and suppose that 



VReuo 
Vim no 



Re uq 
Im no 



and GUi 



and Ui 



VReni 
Vlmtti 



Reni 
Imni / ' 



thus we have 

VRenoVReni + VlmnoVImni = GUq ■ GUj . 
We also use the relation 

|VnoP = \GUo\^. 

We solve problem ()3ip -(j32 p for the real and imaginary parts of i and using our 
notations we obtain the system 



div 



70 



(gU,-2^(^.GU, 



V ^ \GUo\ \\GUo 
+ dw{EoGUi) + diY{EoGUo) + k'^qUi = 0, 
Ui\dn = 0. 

The vector Oq = ^q^"^^ is a unit vector. We can rewrite our system in the form 

div [70 (Id - 200 ® ^o) GUi] + dw{EoGUi) + div(^oGJ7o) + k^qUi = 0, 
or using eigenvectors 



div 



70 ( 6*0 ^0 - ^0 fo ) GUi + dw{EoGUi) + dw{EoGUo) + k^qUi = 0. 
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We suppose that GUi \\ Oq and obtain 

-div(7oVReni) + div(^oVReni) + div(£;oV Reuo) + k'^qReui = 0, (33) 
Reuilao = 0; 

- div (70V Im -ui) + div{EoV Im tti) + div(£'o V Im uq) + k'^q Im ui = 0, (34) 
Imttilao = 0. 

This gives ui. 

Step 4. We calculate 

7 = 70 + ^71 = —}—n(x) -27o(VRenoVReVui + VImMoVImVui)). (35) 
|Vnor 

We set now 70 = 7, and return to the first step to find the corresponding uq and 
repeat the procedure. 



3.2. Algorithm of reconstruction of q for a constant 7 

Step 1. We start from an initial guess qq, and solve the corresponding Dirichlet 
problem for the Helmholtz equation 

7Ano + k'^qo{x)uo = 0, 
uo\dQ = ^• 

Solving the direct problem for tp = arctan y/x ^ obtain uq. 

Step 2. We have seen that our inverse problem is asymptotically approached by 
the direct problem 

|7An + a;2|M^ = 0,in 17, ^^^^ 
\ u = ^p, on do,. 



We compute the difference 



and verify 



eo := 1 — To - Qo (37) 

■"0 



kol < Cprcc, (38) 



where Cprec is our wished order of precision. If condition ()38p holds, we finish our 
algorithm and set q = q^. Otherwise we go to the next step. 
Step 3. We use now the expression 



{qo + 6qi)\uo + 6ui\^ = j{x), 



having the goal to approximate the known function j{x) with the help of the small 
correctors 6ui and 6qi. 
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By expanding the expression, we obtain 

X f, , f^euoReui + lmuolmui\ ,2Wi\'^\ i(^) 

(5gi 1 + 2(5 . — f2 + ^ ] — f2 = 1 — f2 ~ 90 

V V I'^ol / \uo\ J For 

-2(5 . — 5^qo-. — r^. 

For For 

As in Section [3 -H we suppose that (5 <C 1 and that 6 max \qi\ and (5 max |ni| are of 

X X 

the order of 6. Consequently, we consider only terms of order not smaller than 6: 

r Ji^) qo{ReuoReui+lmuolmui) 
^Qi = 1 — - qo - 2(3 , — p5 . 

For For 

To find the corrector ui = 6ui, we expand the following equation 

7A(tio + Sui) + k'^iqo + dqi){uo + 6ui) = 0. 

Considering the terms of order no smaller than 6 and replacing 6qi by the ap- 
proximated formula, we find ui as a solution of the following problem 

'jAui + k -. — -KUi — 2k uQ j — = -eoK no, (39) 

For For 

■Ul\dQ. = 0. 

We solve the problem and obtain ui. 
Step 4. We calculate 

q = qo + 5qi = (j(x) - 2go(Re uq Re ui + Im uq Im ui)) . (40) 

For 

We set now go = Qj and return to the first step to find the corresponding uq and 
repeat the procedure. 



4. Numerical results 



To study the efficiency of this approach, we have tested this method on various 
problems and domains, using the partial differential equation solver FreeFem++ 0. 
We present here one such test. The domain O is a disk of radius 8 centred at the 
origin, which contains three inclusions: a triangle, an L-shaped domain and an el- 
lipse, which represents a convex object, a non-convex object, and an object with 
a smooth boundary respectively. In Figure [T] (a) (respectively (b)) the background 
conductivity (respectively permittivity) is equal to 1 (respectively 3), the conduc- 
tivity (respectively permittivity) takes the value 2.5 (respectively 2) in the triangle, 
1.75 (respectively 1) in the ellipse and 3.05 (respectively 2.55) in the L-shaped do- 
main for the two frequencies /ci = vr • 10^ and /c2 = • 10~^. We purposely choose 
values corresponding to small and large contrast with the background. The initial 
guess in Figure [U (c) (respectively (d)) is equal to 3.5 (respectively 11.5) inside 
the disk of radius 6 centred at the origin, and equal to the supposedly known 
conductivity (permittivity) 1 (respectively 3) near the boundary (outside the disk 
of radius 6). Figure [2] shows the result of the reconstruction when perfect mea- 
sures (with "infinite" precision) are available. For all presented numerical results 
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Figure 1. (a) Distribution of the conductivity 7. (b) Distribution of the permittivity q. (c) Initial guess 
for 7. (d) Initial guess for q. 




< J 



(a) 



(b) 



(c) 




(d) 



Figure 2. Reconstruction test with a "perfect" mesh, (a) Collected data J for the reconstruction of 7. 
(b) Collected data j for the reconstruction of q. (c) Reconstructed conductivity 7. (d) Reconstructed 
permittivity q. 



we use as boundary potential ip = e 



i arctan(z/?/) 



. Figure [2] (a) and (b) represents the 



collected data, J{x) and j{x). For known values of the contrast, we remark that 
through we can 'see' the structure of the permittivity. On Figures [2] (c) and (d), the 
reconstructed conductivity and permittivity are represented: they perfectly match 
the target. Figure H] (a) (respectively H] (b)) presents different errors as functions 
of the iteration for 7 (respectively q). 




(a) 



(b) 




(c) 



(d) 



Figure 3. "Imperfect" meshes for: (a) 50 boundary points; (b) 100 boundary points; (c) 200 boundary 
points; and (d) 400 boundary points. 

We have also considered imperfect data. We ran the reconstruction algorithm 
with the same conditions, but now assume that the data was measured at the 
nodes of a regular mesh on the disk, with 50, 100, 200 and 400 boundary points 
(see meshes on Figure [3]). Figure [5] shows the obtained reconstructions, which still 
perfectly match the target. The convergence result for different number of boundary 
points is given on Figure[6]for the errors — qWl^ and ||J/|Vnp — 7||l^- We 

can observe that the convergence is exponential and that it is even more faster for 
meshes of 50 and 100 boundary points than for meshes of 200 or 400 boundary 
points. 

This better convergence for meshes with 50 and 100 boundary points can be 
illustrated by the following example. For all types of mesh we can perfectly recon- 
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Figure 5. Reconstruction tests for different "imperfect" meshes: (a) 7 and (b) q using a regular mesh with 
50 boundary points, (c) 7 and (d) q using a regular mesh with 100 boundary points, (e) 7 and (f) q using 
a regular mesh with 200 boundary points and (g) 7 and (h) q using a regular mesh with 400 boundary 
points. 

struct 7 and q by the perturbative method if one of the chosen frequency (for the 
reconstruction of 7) is big enough and the second frequency (for the reconstruction 
oi q) is small enough. In the previous examples, the frequencies were chosen equal 
to /ci = vr X 10^ and k2 = t^x 10^^. During our numeric simulations, we have noticed 
that the smaller \ki — /C2I becomes, less efficient the convergence. More precisely, 
the algorithm does not converge for — /E2I < 10 for the case of meshes with 200 
and 400 boundary points (see Figure ([7])). 

Let us analyse the explanation of these results. We notice that there are two 
necessary conditions to be satisfied to ensure the convergence of the algorithm by 
perturbations: 

(1) Using the approximation 7(x) = ^■^^^^-^^2 and q{x) = ji^\2, we need to 
ensure that there exist 5i > and ^2 > such that for each iteration step 
for / = 1,2,..., I Vn^ (x)]"^ > 5i and (x)|^ > 62 (where ki are the chosen 
frequencies). In other words, we need that the sequences {|Vtt'j,.(a;)p : I = 
1,2...} and {\u[^{x)\'^ : / = 1,2...} have some uniform positive lower 
bound. 

(2) The corrector functions to update the initial guess for 7 and q should be 
small enough {\ui \ = 6\ui\ ^ 1) and for I 00, \u\\ should tends to 0. 
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10 

Number of iterations 



Figure 6. Convergence results. Errors Hj/ltip 
number of boundary points on q and 7. 



and II J/|Vup — 7||l^ for meshes with different 
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Figure 7. Dependence of the errors || Jfej/|Vjip — 7||loc s-^d llife2/l^'^P ~ ^lUcxj values of fci, fc2- 

(a) Case of 50 points on the boundary, (b) Case of 100 points on the boundary, (c) Case of 200 points on 
the boundary, (d) Case of 400 points on the boundary. 

Indeed, if the first condition does not hold, we have a division by zero and the 
algorithm has no any sense. In the second condition, the smallness of the correctors 
functions u\ is the basic assumption for deriving the approximate systems (I33p -(j34p 
and (j39p which avoid all the terms of the second order on S. If \u\\ is not small 
enough, we cannot do it any more and the solutions of systems ([33]) - p^ and ([39]) 
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Number of iterations 



(c) 

Figure 8. Plot of max P versus the number of iterations for different values of fci = tt X 10"*, where iii 
is the corrector in the reconstruction of 7. (a) Case of 50 points on the boundary, (b) Case of 100 points 
on the boundary, (c) Case of 200 points on the boundary. 

have no any sense. Moreover, the algorithm converges if and only if l-u'^^l — t- for 
/ — oo. 

Figure [8] shows the decay behaviour of the upper bound of l^xP for the corrector 
Ui from the conductivity update algorithm (see system ([33|) - ([M|) ) for different 
frequencies and meshes. We observe that we have a good convergence corresponding 
to the logarithmic decay of \ui\'^ for all frequencies and meshes with 50 and 100 
boundary points, but we have a divergence result corresponding to the non-decay 
of |ttiP for the frequency ki = lOvr and for the mesh with 200 boundary points. 
The corrector function u[ for the reconstruction of q in our numerical tests for 
/c2 = vr X 10"™, m = 1, 2, 3, 5, 7, is equal to zero. This means that at each iteration 
step we update qo by |:^^r^^- 

To understand why for ki = lOvr and k2 = O.lvr the algorithm diverges for a mesh 
of 200 boundary points and converges for a mesh of 50 or 100 boundary points, let 
us verify the first condition of the convergence. Figure [9] (respectively Figure [TOj) 
shows the lower bounds of |Vn^.(x)p and |tifc.(2;)p for different ki (respectively k2) 
and for different meshes. We notice that we have for all cases, except the case for 
ki = lOvr, k2 = O.lvr and for the mesh with 200 boundary points, that the sequences 
{min iVni (x)\'^ : / = 1, 2 . . .} and {min lui (x)\'^ : / = 1, 2 . . .} converge for / — ?> oo 

a; ' a; » 

to a positive constant. Therefore, we see that for ki = lOvr and /c2 = O.Itt we 
obtain a divergence of the quantities y^^'^^^^p and j:;;r^^- The divergence does 

not take place for a small number of boundary points because of a lower order 
of the precision. For example, we notice that with the growth of the number of 
boundary points {i.e. with the growth of the precision) the limits of the sequences 



April 21, 2010 



0:14 Applicable Analysis Helm 



18 



Taylor & Francis and I. T. Consultant 




Number of iterations 

(a) 




Number of iterations 

(b) 



C 

6 2 







■ ■ - ■ - m 


= 1 


f , , m 


= 2 


* — ■ — m 


= 3 


m 


= 5 


' — »— m 


= 7 










yvw — 





Number of iterations 

(c) 



. B - m — 1 

-•- ■ m — 2 
— ■ — m — 3 
m — 5 
— « — m — 7 



Number of iterations 

(d) 




Number of iterations 

(e) 




Number of iterations 

(f) 



Figure 9. Plot of min|Vit|^ and min|?i|^ versus the number of iterations for different values of fci = 
TT X 10™, where u is the numeric solution of the Helmholtz problem, (a) and (b) Case of 50 points on the 
boundary, (c) and (d) Case of 100 points on the boundary, (e) and (f) Case of 200 points on the boundary. 



{mill I Vn^ 



/ = 1, 2 . . .} and {min {x)\'^ '■ I = 1,2...} becomes more and 



more smaller, as illustrated on Figure [TT] for min \ui P. In the case the mesh of 400 

X 1 

boundary points, the divergence stops the numeric test by the error of the division 
by zero. 



Appendix A. Proof of Lemma II. II 



Let D{x) 
and /(x) 



F{u)f{xa) + G{u){bx — 1), where a 

Using the linearity of the second term, and by introducing 



- with b 

7 



I are unknown 



x+l 
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Figure 10. Plot of min|Vup and min|itp versus the number of iterations for different values of k2 = 
TT X 10~™, where u is the numeric solution of the Helmholtz problem, (a) and (b) Case of 50 points on the 
boundary, (c) and (d) Case of 100 points on the boundary, (c) and (f) Case of 200 points on the boundary. 

N{x) = F{u)f{xa) — D{x), we see that 

^^^^ ^ N{x2)-N{xi) ^ ^ X2N{xi) -XiNix2) 

X2 - Xl X2 - Xl 

_ N{xi)ix2 -X) + N{X2){X - Xl) 
X2 - Xl 

By returning to D, and by introducing the function 

D{xi){x2 - X) + D{X2){X - Xl) 



d{xi,X2,x) = D{x) 



X2 - Xl 
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Figure 11. Plot of min |«|^ versus the number of iterations for different values of fci = vr X 10™, m = 1,5, 
where u is the numeric solution of the Helmholtz problem. 



we have 

d{i,X2,x) = F{u) 
which is also 
d{xi,X2^x) = F{u) 4a^ 

Let us define 

Q{xi,X2,x:i,a) = 4a^ 

We have obtained 



f{xa) 



f{axi){x2 - xi) + f{ax2){x - xi) 



X2 - Xi 



X[X — Xl — X2) + X1X2 



a^xxiX2 + a?{x{xi + X2) + X1X2) + a{xi + X2 + x) + 1 



X3(X3 - Xi - X2) + X1X2 



0^x3X1X2 + a2(x3(xi + X2) + X1X2) + a(xi + X2 + X3) + 1 



d{xi,Xj,Xk) = F{u)Q{xi,Xj,Xk,a). 

Note that Q{xi,Xj,Xk,a) = Q{xj,Xi,xi.,a), but other permutation do not in gen- 
eral yield the same values. 

As a consequence, from n distinct measurements, we obtain 3C^ identities, that 
is, 3C^ formulas of the form 



Fiu) 



Q{xi, Xj , Xfc, fl) 



d{xi , Xj , x/j) 



(Al) 



The value of a can thus be deducted by intersection. 

Note that Q, as a function of a, has only two roots equal to zero (for X3(x3 — 
xi — X2) + X1X2 / 0). By an appropriate choice of Xj, Xj, x^, we can set a G (0, 00). 

We see that the equation becomes 

^/ \ d(^Xi,Xj,Xjf) . / / / X 

Iq/yXi, Xj , X};, a) — -— -. -. Y—i^[X^,Xj,Xj^,Cl). 

d(x'i,x;.,x'j 

Provided that the function a 1-^ Qi^y^j,^k,a) ^jjjgg^jve, a is determined uniquely. 

^ ^x j ^x ^ ^dj 

By using relation ()Aip . this defines F, and therefore N and finally G. 
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Consequently, to determine a and b, it is sufficient to choose four different points 
xi, X2, X3 and X4 to obtain a bijective function on (0, 00) of the form 
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